!> author: 左志华
!> date: 2022-06-07
!>
!> Variable-length array structure <br>
!> linear array vector data structure <br>
!> 可变长度数组结构体定义:
!> 线性数组向量数据结构
!> @todo 视情况，待算法成型，考虑优化，即暂不优化
module sph_list_type

    !> list <br>
    !> 列表
    type list_t
        integer, allocatable :: items(:)    !! list items <br>
                                            !! 表元素
        integer :: log                      !! list length <br>
                                            !! 表元素个数
    contains
        procedure, private :: expand => list_t_expand
    end type list_t

contains

    !> init length <br>
    !> 初始化长度
    elemental subroutine list_t_init(self, capacity)
        type(list_t), intent(inout) :: self     !! list <br>
                                                !! 表
        integer, intent(in) :: capacity         !! list capacity <br>
                                                !! 向量长度

        if (.not. allocated(self%items)) then
            allocate (self%items(capacity))
        else if (size(self%items) /= capacity) then
            deallocate (self%items)
            allocate (self%items(capacity))
        end if

        self%log = 0

    end subroutine list_t_init

    !> Push item
    !> 推送值
    pure subroutine list_t_push(self, value)
        type(list_t), intent(inout) :: self     !! list
                                                !! 表
        integer, intent(in) :: value            !! value
                                                !! 值

        self%log = self%log + 1
        if (self%log > size(self%items)) then
            call self%expand()
        end if
        self%items(self%log) = value

    end subroutine list_t_push

    !> Merge list
    !> 合并列表
    pure subroutine list_t_merge(self, values)
        type(list_t), intent(inout) :: self     !! list
                                                !! 表
        type(list_t), intent(in) :: values      !! values
                                                !! 值

        do while (self%log + values%log > size(self%items))
            call self%expand()
        end do

        self%items(self%log + 1:self%log + values%log) = values%items(:values%log)
        self%log = self%log + values%log

    end subroutine list_t_merge

    !> Expand list size
    !> 拓展空间
    !> @note 空间低于 128 拓展到 128；高于 128，拓展到当前数组长度的 2 倍数
    pure subroutine list_t_expand(self)
        class(list_t), intent(inout) :: self    !! list <br>
                                                !! 表
        integer :: n128

        n128 = size(self%items)/128
        if (n128 == 0) then
            self%items = expand_array(self%items, 128)      ! @tofix vtune 剖析结果：耗时
        else
            self%items = expand_array(self%items, 2*size(self%items))
        end if

    contains

        !> Expand array
        !> 拓展数组
        pure function expand_array(in, capacity) result(out)
            integer, intent(in) :: in(:)
            integer, intent(in) :: capacity
            integer, allocatable :: out(:)

            allocate (out(capacity))
            out(1:size(in)) = in

        end function expand_array

    end subroutine list_t_expand

end module sph_list_type

!> Twin list structure <br>
!> 孪生数组向量数据结构
!> @note 用于存储粒子对
module twin_sph_list_type

    use sph_kinds, only: rk

    !> Twin list <br>
    !> 孪生列表
    type twin_list_t
        integer, allocatable :: items(:, :)     !! list items <br>
                                                !! 表元素
        real(rk), allocatable :: w(:)           !! Smooth kernel value <br>
                                                !! 光滑核函数值
        real(rk), allocatable :: dwdx(:, :)     !! Smooth kernel derivative <br>
                                                !! 光滑核函数导数值
        integer :: log                          !! list length <br>
                                                !! 表元素个数
    contains
        procedure, private :: expand => twin_list_t_expand
    end type twin_list_t

contains

    !> init length <br>
    !> 初始化长度
    pure subroutine twin_list_t_init(self, capacity, save_capacity)
        type(twin_list_t), intent(inout) :: self    !! list <br>
                                                    !! 表
        integer, intent(in) :: capacity             !! list capacity <br>
                                                    !! 向量长度
        logical, intent(in) :: save_capacity        !! save capacity <br>
                                                    !! 是否保存扩展后的长度

        if (.not. allocated(self%items)) then
            allocate (self%items(2, capacity), self%w(capacity), self%dwdx(2, capacity))
        else if (.not. save_capacity .and. size(self%items, 2) /= capacity) then
            deallocate (self%items, self%w, self%dwdx)
            allocate (self%items(2, capacity), self%w(capacity), self%dwdx(2, capacity))
        end if

        self%log = 0

    end subroutine twin_list_t_init

    !> Push item <br>
    !> 推送值
    pure subroutine twin_list_t_push(self, items, w, dwdx)
        type(twin_list_t), intent(inout) :: self    !! list <br>
                                                    !! 表
        integer, intent(in) :: items(2)             !! values <br>
                                                    !! 值
        real(rk), intent(in) :: w                   !! Smooth kernel value <br>
                                                    !! 光滑核函数值
        real(rk), intent(in) :: dwdx(2)             !! Smooth kernel derivative <br>
                                                    !! 光滑核函数导数值

        self%log = self%log + 1
        if (self%log > size(self%items, 2)) then
            call self%expand()
        end if
        self%items(:, self%log) = items
        self%w(self%log) = w
        self%dwdx(:, self%log) = dwdx

    end subroutine twin_list_t_push

    !> Expand list size <br>
    !> 拓展空间
    pure subroutine twin_list_t_expand(self)
        class(twin_list_t), intent(inout) :: self   !! list <br>
                                                    !! 表

        associate (len => 2*size(self%items, 2))
            self%items = expand_matrix_int(self%items, len)
            self%w = expand_vector_real(self%w, len)
            self%dwdx = expand_matrix_real(self%dwdx, len)
        end associate

    end subroutine twin_list_t_expand

    !> Expand matrix size <br>
    !> 拓展数组
    pure function expand_matrix_int(in, capacity) result(out)
        integer, intent(in) :: in(:, :)
        integer, intent(in) :: capacity
        integer, allocatable :: out(:, :)

        allocate (out(2, capacity))
        out(:, 1:size(in, 2)) = in

    end function expand_matrix_int

    !> Expand matrix size <br>
    !> 拓展数组
    pure function expand_matrix_real(in, capacity) result(out)
        real(rk), intent(in) :: in(:, :)
        integer, intent(in) :: capacity
        real(rk), allocatable :: out(:, :)

        allocate (out(2, capacity))
        out(:, 1:size(in, 2)) = in

    end function expand_matrix_real

    !> Expand vector size <br>
    !> 拓展向量
    pure function expand_vector_real(in, capacity) result(out)
        real(rk), intent(in) :: in(:)
        integer, intent(in) :: capacity
        real(rk), allocatable :: out(:)

        allocate (out(capacity))
        out(1:size(in)) = in

    end function expand_vector_real

end module twin_sph_list_type

!> author: 左志华
!> date: 2022-07-16
!>
!> Background grid list (2-dimensional) <br>
!> 背景网格链表法 (2 维)，为所有粒子（包括实、虚粒子）进行背景网格匹配（experimental code）
module sph_linked_list_nnps

    use sph_kinds, only: rk
    use sph_math, only: get_distance, const
    use sph_list_type, only: list_t, list_init => list_t_init, list_push => list_t_push, merge => list_t_merge
    use twin_sph_list_type, only: twin_list_t, twin_init => twin_list_t_init, twin_push => twin_list_t_push
    use sph_smoothed_kernel_function, only: skf
    implicit none
    private

    public :: query, build_grid_from_loc, estimated_capacity
    integer :: estimated_capacity   !! 预估粒子对数目

contains

    !> Build background grid (build storage structure) <br>
    !> 建立背景网格 (构建存储结构)
    pure subroutine setup_grid(scale_hsml, min, max, estimated_capacity, grid)
        real(rk), intent(in) :: scale_hsml                            !! smoothing length <br>
                                                                !! 光滑长度
        real(rk), intent(in) :: min(2)                          !! minimum value <br>
                                                                !! 最小值
        real(rk), intent(in) :: max(2)                          !! maximum value <br>
                                                                !! 最大值
        integer, intent(in) :: estimated_capacity               !! estimated capacity <br>
                                                                !! 预估表存储元素量
        type(list_t), intent(out), allocatable :: grid(:, :)    !! grid <br>
                                                                !! 网格
        integer :: i, j

        i = int((max(1) - min(1))/scale_hsml) + 1
        j = int((max(2) - min(2))/scale_hsml) + 1

        allocate (grid(i, j))
        call list_init(grid, estimated_capacity)

    end subroutine setup_grid

    !> Deploy particle to background grid (mapping) <br>
    !> 部署粒子至背景网格 (映射粒子，mapping)
    pure subroutine mapping_to_grid(i, loc, min, scale_hsml, grid)
        integer, intent(in) :: i                    !! particle index <br>
                                                    !! 粒子编号
        real(rk), intent(in) :: loc(2)              !! particle location <br>
                                                    !! 粒子坐标
        real(rk), intent(in) :: min(2)              !! minimum value <br>
                                                    !! 最小值
        real(rk), intent(in) :: scale_hsml                !! smoothing length <br>
                                                    !! 光滑长度
        type(list_t), intent(inout) :: grid(:, :)   !! grid <br>
                                                    !! 背景网格
        integer :: index(2)

        index(:) = int((loc(:) - min(:))/scale_hsml) + 1

        call list_push(grid(index(1), index(2)), i)

    end subroutine mapping_to_grid

    !> Construct background grid from particle location <br>
    !> 从粒子位置构建背景网格
    pure subroutine build_grid_from_loc(loc, scale_hsml, estimated_capacity, grid)
        real(rk), intent(in) :: loc(:, :)                       !! particle location <br>
                                                                !! 粒子坐标
        real(rk), intent(in) :: scale_hsml                            !! length <br>
                                                                !! 边长
        integer, intent(in) :: estimated_capacity               !! estimated capacity <br>
                                                                !! 单个网格内的预估表存储元素量
        type(list_t), intent(out), allocatable :: grid(:, :)    !! grid <br>
                                                                !! 背景网格
        integer :: i

        associate (min => minval(loc, dim=2) - const(1))    ! const(1) 是小量，用于扩宽区域
            call setup_grid(scale_hsml, min, maxval(loc, dim=2), estimated_capacity, grid)
            do concurrent(i=1:size(loc, 2))
                call mapping_to_grid(i, loc(:, i), min, scale_hsml, grid)
            end do
        end associate

    end subroutine build_grid_from_loc

    !> Query twin particle list <br>
    !> 搜索 (query)
    !> @note L 型查找方式，借鉴 SPHysics 的方案；
    !> 使用查询功能，需要先指定核函数指针内容
    !> 不能查找 1*1 和 2*2 的网格，必须是 3*3 及以上的网格！
    !> @todo 优化：采用数组列优先原则，理论上能稍微提升性能
    pure subroutine query(loc, grid, scale, estimated_capacity, twin_list)
        real(rk), intent(in) :: loc(:, :)               !! particle location <br>
                                                        !! 粒子坐标
        type(list_t), intent(in) :: grid(:, :)          !! grid <br>
                                                        !! 网格
        integer, intent(in) :: scale                    !! scale <br>
                                                        !! 尺度
        integer, intent(in) :: estimated_capacity       !! estimated capacity <br>
                                                        !! 预计交互的粒子对数量
        type(twin_list_t), intent(inout) :: twin_list   !! twin list <br>
                                                        !! 粒子对
        type(list_t) :: found
        integer :: i, j, m, n

        call twin_init(twin_list, estimated_capacity, .true.)
        m = size(grid, 1); n = size(grid, 2)
        ! 仅处理中间区域 (\[i \in [2, size(grid, 1)-1],  j \in [2, size(grid, 2)]\])
        ! |-------------|
        ! | |---------| |
        ! | |         | |
        ! | |         | |
        ! | =========== |
        do j = 2, n
            do i = 2, m - 1
                call list_init(found, 64)
                if (grid(i, j)%log == 0) cycle  ! 网格内无粒子
                ! 查询项的 L 型网格
                call merge(found, grid(i - 1, j - 1))   ! x x x
                call merge(found, grid(i, j - 1))       ! _ o x
                call merge(found, grid(i + 1, j - 1))   ! _ _ _
                call merge(found, grid(i + 1, j))

                call find_nearby_particles(grid(i, j), found, twin_list)
            end do
        end do

        ! i = [1], j = 2, n
        i = 1
        do j = 2, n            ! (1, 2:n) 单元格 （左边）
            call list_init(found, 64)
            if (grid(i, j)%log == 0) cycle
            call merge(found, grid(i, j - 1))       ! x x
            call merge(found, grid(i + 1, j - 1))   ! o x
            call merge(found, grid(i + 1, j))       ! _ _

            call find_nearby_particles(grid(i, j), found, twin_list)
        end do

        i = m
        j = 1
        if (grid(i, j)%log /= 0) then   ! (m, 1) 单元格 （点），无相邻 L 型网格
            call list_init(found, 0)
            call find_nearby_particles(grid(i, j), found, twin_list)
        end if

        do j = 2, n            ! (m, 2:n) 单元格 （右边）
            call list_init(found, 64)
            if (grid(i, j)%log == 0) cycle
            call merge(found, grid(i, j - 1))       ! x x
            call merge(found, grid(i - 1, j - 1))   ! _ o

            call find_nearby_particles(grid(i, j), found, twin_list)
        end do

        ! i = 1, m - 1, j = [1]
        j = 1
        do i = 1, m - 1        ! (1:m-1, 1) 单元格 （上边）
            call list_init(found, 64)
            if (grid(i, j)%log == 0) cycle
            call merge(found, grid(i + 1, j))   ! _ o x

            call find_nearby_particles(grid(i, j), found, twin_list)
        end do

    contains
        !> Find nearby particles <br>
        !> 查找邻近粒子 (find_nearby_particles)
        !> - 相同网格内；
        !> - L 型邻网格。
        !> @note 简化代码
        pure subroutine find_nearby_particles(grid, found, twin_list)
            type(list_t), intent(in) :: grid                !! grid <br>
                                                            !! 中心网格粒子列表
            type(list_t), intent(in) :: found               !! nearby particles <br>
                                                            !! L 型网格粒子列表
            type(twin_list_t), intent(inout) :: twin_list   !! twin list <br>
                                                            !! 粒子对
            integer :: ii, ij
            real(rk) :: dx(2), r, w, dwdx(2)
            do ii = 1, grid%log
                do ij = ii + 1, grid%log
                    call get_distance(loc(:, grid%items(ii)), loc(:, grid%items(ij)), dx, r)
                    if (r < skf%scale_hsml) then
                        call skf%kernel(r, dx, w, dwdx)
                        call twin_push(twin_list, [grid%items(ii), grid%items(ij)], w, dwdx)
                    end if
                end do
                do ij = 1, found%log
                    call get_distance(loc(:, grid%items(ii)), loc(:, found%items(ij)), dx, r)
                    if (r < skf%scale_hsml) then
                        call skf%kernel(r, dx, w, dwdx)
                        call twin_push(twin_list, [grid%items(ii), found%items(ij)], w, dwdx)
                    end if
                end do
            end do
        end subroutine find_nearby_particles

    end subroutine query

end module sph_linked_list_nnps
